Índice

Professor: Marcio Argollo

Contato: Sala A1-02, tel. 2629-5773

Referências básicas:

• Manuais de Fortran 77:
  • http://www.idris.fr/data/cours/lang/fortran/f90/F77.html
  • http://www.stanford.edu/class/me200c/tutorial_77/
• Manuais de Fortran 90
  • http://www.cenapad.unicamp.br/servicos/treinamentos/apostilas/apostila_fortran90.pdf (baixar aqui)
  • http://www.cisl.ucar.edu/tcg/consweb/Fortran90/F90Tutorial/tutorial.html
• Fortran 90 para programadores de Fortran 77
  • http://www.nsc.liu.se/~boein/f77to90/
• Manuais de C/C++
  • http://www.cenapad.unicamp.br/servicos/treinamentos/apostilas/apostila_C.pdf (baixar aqui)
  • http://imada.sdu.dk/~svalle/courses/dm14-2005/mirror/c/
  • http://www.cs.cf.ac.uk/Dave/C/CE.html
  • http://www.ic.unicamp.br/~cmrubira/aacesta/cpp/cpp15.html
• Livros-texto
  • An Introduction to Computer Simulation Methods: Applications to Physical Systems, H. Gould e J. Tobochnik.

Para uma revisão de comandos básicos em Unix/Linux, utilização do gnuplot e programação em C veja o tutorial apresentado na Semana da Física 2009 "Ferramentas de Computação em Física".

As listas deverão ser entregues em papel e por email, compactadas no formato tar.gz e com o nome listaN_seunome.tar.gz

1. Somas não são comutativas no computador! Calcule as somas e com i sendo uma variável ponto flutuante de precisão simples e coloque em um mesmo gráfico as diferencas vs N e vs N para onde .
2. Escrever um programa que recebe como entrada um número inteiro em representação decimal e retorna sua expansão binária. Para tal, crie uma função que verifica se o i-ésimo bit de uma palavra (variável do programa) vale zero ou um.
3. Considere um número de ponto flutuante de N bits com E bits reservados para o expoente e M=N-E-1 para a mantissa. Represente TODOS os números obtidos nesta representação em um grafico de barras no gnuplot. Para cada número x, armazene um par (x,1) e plote o arquivo de dados com with impulses. Reveja aqui a definição de números de ponto flutuante.
   
4. Imagine o empilhamento de dominós de massa unitária e comprimento 2. O primeiro dominó se encontra com o centro de massa na origem, e a cada passo n deslocamos os dominós ja empilhados de para a direita e adicionamos outro dominó na parte inferior da pilha, também com seu centro de massa na origem. Escreva um programa que diz, para um dado , quantos dominós podem ser empilhados até que o centro de massa ultrapasse a largura do primeiro dominó e a pilha tombe. (Veja uma discussão do problema aqui).

Sendo a derivada de uma função definida por e , o método de Euler (com difereças posteriores) corresponde ao truncamento desta expansão em ordem h. Obtemos assim que a derivada no ponto x é aproximada por .

1. Mostre que o método de Euler com diferenças centradas é exato em , onde h é o passo de discretização da derivada. Nesse método escrevemos
2. Calcule a derivada numérica de no ponto x=3 usando o método de Euler simples com passo de discretização entre e e faça um gráfico mostrando a diferença relativa entre a derivada numérica e o valor exato da derivada para ponto flutuante com precisão simples e dupla. Entre quais valores de h a aproximação é aceitável? (Solução aqui)
3. Mostre que o método de Euler é instável para qualquer passo h na solução do problema de crescimento exponencial
4. Integre a equação diferencial com com o método de Euler simples para diferentes valores de passo h=0.1, h = 0.05, h = 0.01, h = 0.005 e h = 0.001 e compare cada solução em t=1,2,3,4 e 5 com a solução exata . Voce pode encontrar essa solução exata digitando a equação diferencial no site do WolframAlpha
5. Faça um gráfico com a solução para h=0.05 entre t=0 e t=5 e a solução exata.
6. Resolva o problema acima com o método de Euler com diferenças centradas e os mesmos valores de h. Para obter o primeiro passo use o método de Euler simples com passo de discretização dez vezes menor que o usado no resto do intervalo.

• Referências
  • http://mathworld.wolfram.com/NumericalDifferentiation.html
  • http://davinci.if.ufrgs.br/wiki/index.php/M%C3%A9todo_de_Euler

A aproximação para a derivada, como sugerida no método de Euler, introduz erros sistemáticos que acarretam na destruição de certos invariantes da dinâmica original, como a energia total no movimento do pêndulo simples. Para reduzir os efeitos da aproximação numérica à derivada, pode-se recorrer a métodos com aproximações de ordem superior (o método de Euler é de primeira ordem no passo de discretização h).

Ver: http://mathworld.wolfram.com/Runge-KuttaMethod.html
Atenção
Para sistemas de equações diferenciais, ver: http://www.nsc.liu.se/~boein/f77to90/rk.html

1. Resolva numericamente o problema do pêndulo físico utilizando o método de Runge-Kutta de quarta ordem com passo e ilustre o espaço de fase do sistema. Considere metros, onde é o comprimento da haste de massa desprezível à qual se prende uma massa m. Use a condição inicial (posição vertical para baixo) e velocidades iniciais e , onde é a velocidade crítica que separa os movimentos de oscilação e rotação (Ver o artigo "Comportamento crítico no pêndulo simples" do Paulo Murilo para um estudo desta transição). Considere a evolução por dois períodos do movimento. Ilustre, no mesmo gráfico, a solução do problema na aproximação de pequenos ângulos com e .
2. Faça um gráfico contendo os primeiros 30 segundos da evolução temporal dos valores analítico e numérico da energia total do pêndulo. No caso numérico, utilize o método de Euler e o método de Runge-Kutta de segunda e quarta ordem com passos e .

Ver Reversible multiple time scale molecular dynamics (pdf aqui)

Estes métodos procuram criar esquemas de discretização da evolução Hamiltoniana de um sistema conservativo, de modo que o volume de pontos representativos no espaço de fase seja preservado durante a evolução temporal. Estes esquemas envolvem fatorização do operador Liouvilleano.

1. Use o método de Euler-Cromer e o método de Euler com passo e para resolver numericamente o problema do pêndulo simples na aproximação de pequenas oscilações. Ulitize a condição inicial e . Faça um gráfico contendo a evolução temporal da posição angular por aproximadamente 10 períodos do movimento para a solução analítica e numérica com os métodos de Euler e Euler-Cromer. Faça outro gráfico ilustrando a evolução temporal da energia no mesmo alcance temporal para as 3 soluções.

Nesta sequência de problemas utilize o método de Runge-Kutta de quarta ordem com passo

Considere o pêndulo simples em um meio viscoso, que oferece resistência ao seu movimento com uma força proporcional à sua velocidade . A equação de movimento pode ser escrita como .

1. Faça o gráfico da posição como função do tempo para e com posição angular inicial e . Note que para valores grandes de q devemos recuperar o movimento sem dissipação (conservativo) e, portanto, devemos utilizar um método simplético, como Verlet ou Euler-Cromer, para integrar a equação de movimento.

Considere agora uma força periódica impulsionando o pêndulo com força . A equação de movimento pode ser escrita como .

1. Construa o gráfico com e e , respectivamente. Observe a mudança de comportamento ao variarmos a amplitude A da força externa. Simule o pêndulo por 100 ciclos da força externa e despreze os primeiros 10 ciclos (quanto vale o período da força externa neste exemplo?).
2. Construa o espaço de fases para cada uma das situações descritas acima com a evolução do ponto representativo
3. Faça um gráfico “estroboscópico”, chamado seção de Poincaré, imprimindo os pontos no espaço de fase correspondentes a intervalos de tempo múltiplos do período da força externa .
4. Obtenha uma série de pontos onde corresponde ao valor da velocidade angular quando vale zero pela n-ésima vez durante a evolução temporal do pêndulo. Despreze os valores obtidos para tempos menores que 10 ciclos da força externa. Construa agora outra seção de Poincaré em um gráfico para os mesmos parâmetros do problema anterior.

O pêndulo descrito acima deve apresentar comportamentos distintos para diferentes parâmetros. Em particular, para e , obtemos comportamento caótico para alguns valores da amplitude da força externa A. Verifique!

Leia aqui sobre o problema.

Considere o mapa logístico

1. Faça os gráficos (escala lin-log) da evolução temporal de x a partir de um valor inicial arbitrário para e .
2. Fazer “gráficos de escada” para os mesmos valores de parâmetro e condição incial x=0.6. Como critério de parada use a condição .
3. Fazer “gráficos de escada” para o caso e condição inicial x=0.3 e para a segunda iterada f(f(x)) com mesma condição inicial.
4. Construir o diagrama de bifurcação do mapa logístico na região .
5. Obter os expoentes de Lyapunov do mapa logístico na região .

O mapa logístico com é ergódico, o que significa que podemos substituir médias temporais por médias sobre configurações de um ensemble e, nesse caso, onde é a fração de vezes que a órbita se encontra entre x e x+dx.

1. Calcule, subdividindo o intervalo [0,1] em 100 subintervalos, o histograma para as primeiras 100000 iterações do mapa logístico com e encontre . Esse caso particular admite a solução analítica . Compare os resultados mostrando-os em um mesmo gráfico.
2. Obter o diagrama de bifurcação para os mapas

(a)

(b)

ambos na região

1. Voce encontra neste arquivo a sequência dos primeiros 10 mil dígitos de .

• Calcular a frequência P(i) de ocorrência de cada dígito i na sequência.
• Calcular a frequência P(i,j) de ocorrência consecutiva de cada par de dígitos (i,j) na primeira metade da sequência. Utilize essa informação para adivinhar, na segunda metade da sequência, qual o dígito j seguinte a um dado dígito i da sequência. Qual a sua frequência de acertos? Como ela se compara com a probabilidade medida na primeira metade?

1. A função de correlação é uma das ferramentas tradicionais de análise de sinais. Ela mede o grau de correlação entre o valor do sinal medido em um instante e outro sinal medido em um instante de tempo subsequente. Note que, em sinais não estacionários, a função de correlação pode depender do instante inicial da medida . Não trataremos desses casos aqui.
2. Calcule a função de correlação para uma série de 10000 pontos obtidos com a evolução do mapa logístico com (ponto de acumulação, infinitas órbitas periódicas instáveis emergem neste ponto), , (transição caos-ordem em uma janela de período 3) e (ergodicidade comprovada para esse valor de parâmetro).

Para ler sobre automata celulares:
http://mathworld.wolfram.com/ElementaryCellularAutomaton.html
http://classes.yale.edu/fractals/CA/welcome.html

1. Desenvolva um código para ilustrar a evolução da regra 90 utilizando somente operações lógicas e um bit para cada sítio do sistema, implementando um sistema de 32 ou 64 bits (dependendo do tamanho da palavra no computador utilizado). Implemente condições de contorno periódicas e livres e condição inicial com um bit central 1 e todos os outros 0.
2. Desenvolva um código para receber uma dada regra de Wolfram como parâmetro de entrada e retornar a evolução temporal de uma condição inicial por T passos de tempo em um sistema com N sítios. Utilize uma palavra de computador por sítio. Mostre graficamente as evoluções temporais para as regras 30, 60, 90, 102, 126 e 150 com N=T=200. Utilize como condição inicial um único sítio (central) com estado s=1 e todos os outros com s=0 e implemente condições de contorno periódicas (s[N]=s[0]) e livres (s[N]=s[-1]=0);
3. Para a regra 30, escolha um sítio e acompanhe sua evolução temporal 0,1,0,0,1,….
   

• Calcule o número médio de bits 1 ao longo da evolução com N=10000 e T=100, 1000, 10000 e 100000.
  
• Calcule a função de correlação
  
• Agrupe os bits na evolução em blocos de 8 bits. Transforme o bloco em um inteiro com sinal. Calcule o valor médio do número obtido. Calcule a função de correlação para a série gerada com tal procedimento. Faça um gráfico com a evolução temporal dos números de 8 bits. Para lembrar da representação de inteiros com sinal leia:
  

http://www.cs.cornell.edu/~tomf/notes/cps104/twoscomp.html
http://en.wikipedia.org/wiki/Two's_complement

Números gerados com o computador, uma máquina de Turing determinística, não pode gerar aleatoriedade com seus elementos lógicos. Podemos, no máximo, gerar sequências de números que passem por determinados testes estatísticos de aleatoriedade. A estes números chamamos pseudo-aleatórios e sua obtenção depende crucialmente do aumento de entropia gerado pela perda de bits (informação) nas operações lógicas - para os mais interessados: http://www.springerlink.com/content/jn7x3365386phn46/

Existem geradores de números verdadeiramente aleatórios, baseados na estocasticidade de processos físicos como emissão de fótons em processos radiativos ou espalhamento em divisores de feixes (veja aqui um exemplo). Esses geradores são usados para operações críticas, como criptografia e bingo.

Uma boa fonte de informações a respeito de números pseudo-aleatórios é a biblioteca GSL, que possui, dentre diversas outras aplicações científicas como cálculo de autovalores e solução de sistemas de equações lineares, rotinas que geram números pseudo-aleatórios com diversas receitas. http://www.gnu.org/software/gsl/manual/html_node/Random-Number-Generation.html

1. Usando o gerador congruencial x=(x*A)mod(B) gere pontos em um espaço 2D e visualize seu resultado. Utilize A=16807 e 65539 e B=2147483648. Estes números podem ser gerados sem a operação mod com números inteiros de 32 bits com sinal. Visualize seus resultados normalizando o número aleatório para gerar pontos com alcance [-1…1]
2. Usando o mesmo gerador com A=16807 e 65539, crie pontos em um espaço 3D. Visualize seus resultados e, girando a figura no gnuplot, tente visualizar a regularidade nos pontos gerados com 65539. Em 2D, visualize x versus 9x-6y+z onde x,y,z são a trinca de números aleatórios correspondente a um ponto no espaço 3D. Veja como o gerador 65539 é ruim para amostrar pontos no espaço.
3. Gere pares de pontos (x,y) aleatoriamente distribuídos em 2D (com x e y entre -1 e 1). Calcule a fração de números gerados que satisfazem . Essa razão deve ser aproximadamente a razao entre a área de um círculo de raio unitário e um quadrado de lado 2. Esta razão vale . Multiplique sua fração de números por 4 e voce terá calculado pi pelo Método de Monte Carlo. Faça um gráfico mostrando a convergência do valor estimado com o método de Monte Carlo para o valor aproximado de 3.14159265358979323846 (por exemplo, mostre no gráfico a diferença entre o valor aproximado e o valor estimado) utilizando o gerador com A=16807 e compare o resultado obtido com o gerador A=65539. Quantos números aleatórios devemos gerar para se obter de precisão na sua estimativa com A=16807?
4. Ruína do jogador: Considere dois jogadores, o primeiro com n1 moedas e o segundo com n2 moedas em um jogo de cara ou coroa. Cada vez que um jogador perde, ele entrega uma moeda para o outro jogador e o jogo se repete até um dos jogadores perder todo o dinheiro. Realize esse jogo um milhão de vezes e calcule a fração de vezes que cada jogador perdeu todo o dinheiro. Esse valor deve ser interpretado como a probabilidade de cada jogador iir à falência. Compare seu resultado com os valores teóricos para as e . Para ler mais sobre o assunto clique aqui